$\forall a, b \in \mathbb{R}^{+} \cup \{0\}, a \geq b \iff a^2 \geq b^2$
$\implies$ \begin{align*} a & \geq b \\ a \cdot a \geq & b \cdot a \geq b \cdot b \end{align*} Thus, $a^2 \geq b^2$
$\impliedby$
If $b \ne 0$, Prove by contradiction, we have $a^2 \geq b^2$, assume $a \lt b$. Since $a, b \geq 0$, we have $a \cdot a \lt b\cdot a \lt b \cdot b$, which introduced a contradiction.
On the other hand, if $b = 0$, then since $a \geq 0$, $a = 0$, thus $a \geq b \blacksquare$
Thus, we finished the proof
Given that we have proved $L^2$, defined as the space of square-integrable functions, equipped with the corresponding inner product, is a vector space, it follows that the inner product $\langle \cdot, \cdot \rangle$ satisfies:
$$ \langle f, g \rangle = \int_{\Omega} f(x) \overline{g(x)} \, dx $$for all $f, g \in L^2(\Omega)$, where $\Omega$ is the domain of definition and $\overline{g(x)}$ denotes the complex conjugate of $g(x)$, for $\mathbb{R}$, $\overline{g(x)} = g(x)$.
Thus, according to the definition of the inner product, and more generally, a metric, we have:
$$ \langle f, f \rangle \geq 0 $$for all $f \in L^2(\Omega)$, which indicates that the inner product of any function with itself is always greater than or equal to zero.
Thus, to prove $d(f, f_n) \geq d(f, f_{n+1})$, by Lemma 1, it is equivalent to prove the inequality:
$$ d(f, f_n)^2 \geq d(f, f_{n+1})^2 $$Let $g = f_{n + 1}$
We've proved that \begin{align*} d(f,g)^2 &= \langle f-g, f-g \rangle = \int (f-g)(\overline{f-g}) \\ &= \langle f, f \rangle + \langle g, g \rangle - 2\langle f, g \rangle \\ &= \int f^2 + \int g^2 - 2 \int fg \\ &= \|f\|^2 + \int (\sum_{i = 1}^{n} c_i e_i)(\overline{\sum_{j = 1}^{n} c_j e_j}) - 2\sum_{i = 1}^{n} c_i \langle f, e_i \rangle \\ &= \|f\|^2 + \sum_{i = 1}^{n} c_i^2 - 2 \sum_{i = 1}^{n} c_i \langle f, e_i \rangle \end{align*}
where $c_j = \langle f, e_j \rangle$. Thus, $d(f, f_{n+1})^2 = \|f\|^2 + \sum_{i=1}^{n+1} c_i^2 - 2\sum_{i=1}^{n+1} c_i \langle f, e_i \rangle = d(f, f_n)^2 + c_{n+1}^2 - 2c_{n+1}\langle f, e_{n+1} \rangle$, where the last part
\begin{align*} c_{n+1}^2 - 2c_{n+1}\langle f, e_{n+1} \rangle &= \langle f, e_{n+1} \rangle^2 - 2\langle f, e_{n+1} \rangle^2 \\ &= - \langle f, e_{n+1} \rangle^2 \\ &\leq 0 \\ \end{align*},thus proved $d(f, f_n)^2 \geq d(f, f_{n+1})^2$. By Lemma 1, since the metric defined by inner product is greater or equal to zero, we have proved $d(f, f_n) \geq d(f, f_{n+1})$.
from pylab import *
import math
inf_p_deri_f = lambda theta, n: cos(n * sin(theta/2) * pi)
def dft(n):
thetas = arange(0, 2*pi, 2*pi/128)
f_values = [inf_p_deri_f(theta, n) for theta in thetas]
dft = rfft(f_values)
return dft
wavenumber = range(0, 65)
# wavenumber = range(0, 128)
dft_n_2 = dft(2)
dft_n_4 = dft(4)
dft_n_8 = dft(8)
dft_n_16 = dft(16)
modulus_n_2 = absolute(dft_n_2)
modulus_n_4 = absolute(dft_n_4)
modulus_n_8 = absolute(dft_n_8)
modulus_n_16 = absolute(dft_n_16)
plot(wavenumber, modulus_n_2, 'b-',
wavenumber, modulus_n_4, 'g-',
wavenumber, modulus_n_8, 'r-',
wavenumber, modulus_n_16, 'y-', )
xlabel("Wavenumber", fontsize=20)
ylabel("|C\u2096|", fontsize=20)
legend(["n = 2", "n = 4", "n = 8", "n = 16"],
loc="lower left")
yscale('log')
xscale('log')
import plotly.graph_objects as go
import numpy as np
# Define the dft function
def dft(theta, n):
f_values = np.cos(n * np.sin(theta / 2) * np.pi)
dft = rfft(f_values)
modulus = np.abs(dft)
return modulus
# Generate theta and wavenumber
theta = np.arange(0, 2 * np.pi, 2 * np.pi / 128)
# wavenumber = np.fft.rfftfreq(128, 2 * np.pi / 128)
wavenumber = list(range(0, 65))
# Create traces for each n value
traces = []
n_values = [2, 4, 8, 16]
colors = ['blue', 'green', 'red', 'yellow']
for n, color in zip(n_values, colors):
trace = go.Scatter(x=wavenumber, y=dft(theta, n), mode='lines', name=f'n = {n}', line=dict(color=color))
traces.append(trace)
# Create the figure and update layout
fig = go.Figure(data=traces)
fig.update_layout(
title='DFT Modulus vs Wavenumber',
xaxis_title='Wavenumber',
yaxis_title='|Cn|',
xaxis_type='log',
yaxis_type='log',
xaxis=dict(titlefont=dict(size=20)),
yaxis=dict(titlefont=dict(size=20)),
legend=dict(title='Legend', itemsizing='constant'),
width=800,
height=600
)
# Show the figure
fig.show()
len(modulus_n_16)
65
modulus_n_2
array([2.81954443e+01, 3.68486870e+01, 4.04070998e+01, 3.55441159e+01,
9.38180169e+00, 1.29772243e+00, 1.13070890e-01, 6.83823201e-03,
3.05218340e-04, 1.04939818e-05, 2.86900104e-07, 6.39254453e-09,
1.18390847e-10, 1.85090350e-12, 2.78565342e-14, 9.82836620e-16,
4.37244394e-15, 4.16729171e-15, 3.34591022e-15, 2.55444565e-15,
4.79747557e-15, 6.81830523e-15, 1.08347752e-14, 4.47274214e-15,
2.74290167e-15, 8.61106897e-15, 8.13238095e-15, 8.35279207e-15,
4.89820817e-15, 1.32899742e-15, 2.38105646e-15, 3.75985035e-15,
1.98602732e-15, 1.20403124e-15, 2.16328171e-15, 2.93045478e-15,
3.55705088e-15, 4.32034197e-15, 2.81500046e-15, 2.05959260e-15,
4.02305951e-15, 2.97655733e-15, 4.55603185e-15, 1.72896574e-15,
5.34493553e-15, 5.85016270e-15, 3.13757278e-15, 7.11556975e-15,
2.39391990e-15, 2.83977684e-15, 4.40423080e-15, 3.84370056e-15,
3.04394388e-15, 3.56579899e-15, 2.00378954e-15, 2.43296885e-15,
2.84659712e-15, 1.97205799e-15, 4.59217260e-16, 2.45672369e-15,
4.48542468e-15, 3.61462429e-15, 4.63401539e-15, 4.64837533e-15,
1.77635684e-15])
wavenumber[:10]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
dft_n_2
array([ 2.81954443e+01+0.00000000e+00j, -3.68486870e+01+5.66213743e-15j,
4.04070998e+01-1.23011590e-14j, 3.55441159e+01-1.57651669e-14j,
9.38180169e+00-2.27595720e-15j, 1.29772243e+00-4.38191150e-15j,
1.13070890e-01-3.54037173e-15j, 6.83823201e-03+2.67792516e-15j,
3.05218340e-04+7.15714042e-16j, 1.04939818e-05-1.02295491e-15j,
2.86900104e-07+2.54805814e-15j, 6.39254453e-09+1.63660725e-15j,
1.18390847e-10+1.23641649e-15j, 1.85090345e-12-4.58869618e-16j,
2.78556185e-14+2.25865830e-16j, 6.18512498e-16+7.63812878e-16j,
-4.05067994e-15-1.64628612e-15j, -3.02644833e-15-2.86477411e-15j,
3.31657570e-15-4.42086878e-16j, 2.18804739e-15-1.31819619e-15j,
4.68953051e-17-4.79724636e-15j, 6.04415902e-15-3.15553925e-15j,
9.83570743e-15+4.54436069e-15j, -1.64166778e-17+4.47271201e-15j,
-1.65923114e-15+2.18413864e-15j, 5.52311115e-15+6.60649317e-15j,
-2.75906301e-15+7.65004518e-15j, -8.34483111e-15+3.64594377e-16j,
-1.53912527e-15-4.65011147e-15j, -8.97791610e-16-9.79900184e-16j,
1.76780320e-15+1.59508674e-15j, -2.69518282e-15+2.62153851e-15j,
-1.77635684e-15+8.88178420e-16j, -8.57530856e-16+8.45181674e-16j,
-1.76780320e-15-1.24685991e-15j, -2.65492207e-15+1.24054587e-15j,
1.53912527e-15-3.20682154e-15j, 1.23940375e-15-4.13874777e-15j,
-2.79899101e-15-2.99794745e-16j, 4.31327182e-16+2.01392119e-15j,
3.57545006e-15-1.84422470e-15j, -1.44587345e-16+2.97304356e-15j,
1.65076407e-15+4.24645785e-15j, 1.69270768e-15-3.52226117e-16j,
-4.82085431e-15+2.30818099e-15j, -8.55779757e-16+5.78723117e-15j,
1.12431639e-15-2.92921075e-15j, -6.29942508e-15-3.30886332e-15j,
-2.16656899e-15+1.01824914e-15j, 1.60193355e-15-2.34481159e-15j,
-3.22087537e-16-4.39243766e-15j, 2.28082505e-15+3.09384406e-15j,
2.78152272e-15+1.23641649e-15j, -1.11450625e-15-3.38715194e-15j,
-1.23989366e-15-1.57411443e-15j, 1.05439508e-16-2.43068301e-15j,
2.39643917e-15-1.53629246e-15j, 1.97151323e-15-4.63496429e-17j,
2.35922393e-16-3.93981112e-16j, 2.22044605e-15+1.05124243e-15j,
2.66453526e-15+3.60822483e-15j, 3.55271368e-15+6.66133815e-16j,
-3.55271368e-15+2.97528556e-15j, -3.55271368e-15+2.99760217e-15j,
1.77635684e-15+0.00000000e+00j])
len(wavenumber)
65
From the textbook, we have the following Theroem

rfft is $\frac{n}{2} + 1$ with $n$ intepolation points, this is due to the symmetry while intepolation a real functionfft, the coeficient is not the same with the one provided by the textbook.I tried to come up with a function that takes into $\frac{n}{2} + 1$ complex coefficients, but I did not get a good one.
Thus, for simplicity, I will use fft for finding the intepolation function for now
def dft(n):
thetas = arange(0, 2*pi, 2*pi/128)
f_values = [inf_p_deri_f(theta, n) for theta in thetas]
dft = fft(f_values)
return dft
wavenumber = range(0, 65)
# wavenumber = range(0, 128)
dft_n_2 = dft(2)
dft_n_4 = dft(4)
dft_n_8 = dft(8)
dft_n_16 = dft(16)
The following are the original function I tried to come up with the intepolation function. Yet due to the optimization, even with the correct number of points intepolation, the intepolation function's coefs is still different from the one provided by the textbook
import numpy as np
def P(t, a_coefficients, b_coefficients, c, d, n):
summation = 0
for k in range(n):
ak = a_coefficients[k] # k-th coefficient for cosine terms
bk = b_coefficients[k] # k-th coefficient for sine terms
summation += ak * np.cos(2 * np.pi * k * ((t - c) / (d - c))) - bk * np.sin(2 * np.pi * k * ((t - c) / (d - c)))
return (1 / np.sqrt(n)) * summation
def get_trig_inter_func(dft, c, d):
a_coefficients = [i.real for i in dft]
b_coefficients = [i.imag for i in dft]
n = len(dft)
def trig_inter_func(t):
return P(t, a_coefficients, b_coefficients, c, d, n)
return trig_inter_func
# def P(t, a_coefficients, b_coefficients, c, d, n):
# summation = a_coefficients[0] # The zero-frequency term has no sine component and should be added only once.
# for k in range(1, n // 2): # Loop through the positive frequencies.
# # Calculate the cosine and sine terms
# cosine_term = np.cos(2 * np.pi * k * (t - c) / (d - c))
# sine_term = np.sin(2 * np.pi * k * (t - c) / (d - c))
# summation += a_coefficients[k] * cosine_term
# summation += b_coefficients[k] * sine_term
# summation += a_coefficients[-1] * np.cos(np.pi * (t - c) / (d - c)) # Nyquist frequency term
# return (1 / np.sqrt(n)) * summation
# def get_trig_inter_func(dft, c, d, n):
# # Handle the symmetric part of the DFT coefficients for a real-valued signal.
# a_coefficients = [dft[0].real] + [coeff.real for coeff in dft[1:n//2]] + [dft[n//2].real]
# b_coefficients = [0] + [coeff.imag for coeff in dft[1:n//2]] + [0] # Imaginary part for Nyquist frequency is 0.
# def trig_inter_func(t):
# return P(t, a_coefficients, b_coefficients, c, d, n)
# return trig_inter_func
I did some research and find the following optimazation used for fft, mainly from numpy's documentation.
Thus, we need to do the following:
The following is the corrected function
In mathmatic form: Let $n$ be the total number of Discrete Fourier Transform (DFT) results or a given $t_{\text{num}}$ if provided. For a complex number $z$, $\Re(z)$ represents the real part, and $\Im(z)$ represents the imaginary part.
The coefficients are defined as:
If $n$ is even, then:
If $n$ is odd, the last term is included in the series and is doubled:
The trigonometric interpolation function $f(t)$ at a point $t$ is given by:
If $n$ is even, include the Nyquist term:
def get_trig_interpolation_function(dft_results, c, d, t_num: int = None):
if t_num is None:
n = len(dft_results)
else:
n = t_num
a0 = dft_results[0].real / n # The constant term (DC component) should not be doubled
an = (dft_results[1:n // 2].real * 2 / n).tolist() # The cosine coefficients
bn = (-dft_results[1:n // 2].imag * 2 / n).tolist() # The sine coefficients
# If n is even, the Nyquist frequency term should be handled separately and should not be doubled
if n % 2 == 0:
an.append(dft_results[n // 2].real / n)
else: # If n is odd, the last term is still a part of the series and should be doubled
an.append(dft_results[n // 2].real * 2 / n)
# Function to calculate the trigonometric interpolation at a given point t
def trig_interpolation_function(t):
result = a0
for k in range(1, n // 2):
result += an[k - 1] * np.cos(2 * np.pi * k * t / (d - c))
result += bn[k - 1] * np.sin(2 * np.pi * k * t / (d - c))
if n % 2 == 0: # Add the Nyquist term
result += an[-1] * np.cos(np.pi * (n // 2) * t / (d - c))
return result
return trig_interpolation_function
x = np.linspace(0, 2 * np.pi, 1000)
thetas = arange(0, 2*np.pi, 2*np.pi/128)
fig = go.Figure()
y = [inf_p_deri_f(i, 16) for i in x]
trig_inte_f = get_trig_interpolation_function(dft_n_16, 0, 2*pi, 128)
s = [trig_inte_f(i) for i in x]
inte_ps = [inf_p_deri_f(i, 16) for i in thetas]
# Add a scatter plot to the figure
fig.add_trace(go.Scatter(x=x, y=y, mode='lines', name='truth'))
fig.add_trace(go.Scatter(x=x, y=s, mode='lines', name='inte'))
fig.add_trace(go.Scatter(x=thetas, y=inte_ps, mode='markers', name='inte points'))
# Update layout for a better view
fig.update_layout(
title='Plot',
xaxis_title='x',
yaxis_title='y',
autosize=False,
width=800,
height=600
)
# Show the plot
fig.show()
x = np.linspace(0, 2 * np.pi, 1000)
thetas = arange(0, 2*np.pi, 2*np.pi/128)
fig = go.Figure()
y = [inf_p_deri_f(i, 16) for i in x]
trig_inte_f = get_trig_interpolation_function(dft_n_16, 0, 2*pi, 128)
s = [trig_inte_f(i) for i in x]
inte_ps = [inf_p_deri_f(i, 16) for i in thetas]
# Add a scatter plot to the figure
fig.add_trace(go.Scatter(x=x, y=y, mode='lines', name='truth'))
fig.add_trace(go.Scatter(x=x, y=s, mode='lines', name='inte'))
fig.add_trace(go.Scatter(x=thetas, y=inte_ps, mode='markers', name='inte points'))
sliders = []
slider_a = dict(
active=0,
currentvalue={"prefix": "Terms used: "},
pad={"t": 50},
steps=[{
'method': 'restyle', # Using 'restyle' to update existing traces
'label': str(a),
'args': [{
'y': [
y, # Reapply existing data for 'truth' trace
list(get_trig_interpolation_function(dft_n_16, 0, 2*pi, a)(i) for i in x), # Update data for 'inte' trace
inte_ps # Reapply existing data for 'inte points' trace
]
},
[0, 1, 2] # Indices of the traces to update: 0 for 'truth', 1 for 'inte', 2 for 'inte points'
]
} for a in range(1, 129)]
)
sliders.append(slider_a)
# Update layout for a better view
fig.update_layout(
sliders=sliders,
title='Plot',
xaxis_title='x',
yaxis_title='y',
autosize=False,
width=800,
height=600,
xaxis=dict(range=[0 - 0.1, 2*pi + 0.1]), # Set fixed range for x-axis
yaxis=dict(range=[-1 - 0.1, 1 + 0.1])
)
# Show the plot
fig.show()
Assume the "wavenumber" refers to the frequency domain representation of a given function. Specifically, it is related to the spatial frequency of a wave, which corresponds to the number of wave cycles per unit distance.
Assume the "modulus of the Discrete Fourier Transform of $f$" refers to the magnitude of the complex numbers that result from applying the DFT to $f$
Let's take a closer look at the equation in the thereom:
The graphs shown ($C_k$ and $k$) illustrate two primary relationships: the magnitude of the Fourier coefficients $C_k$ as a function of the wavenumber $k$, and the critical wavenumber $k_{cut}$ as a function of the parameter $n$.
In the context of Fourier Transformation, the Fourier coefficients $C_k$ can be seen as coordinates in a higher-dimensional space, where each dimension corresponds to a sine or cosine basis function.
As $k$ increases, we are effectively examining higher frequency components of the original function. If the coefficients $C_k$ decay with increasing $k$, it suggests that high-frequency components have less influence on the shape of the original function — in other words, the function is smoother and contains fewer rapid oscillations.
This is akin to a weighted sum of the basis functions, where the weight is given by the magnitude of $C_k$. For lower values of $k$, the weights are larger, signifying that these frequencies are more significant in reconstructing the original function. As $k$ increases (and the corresponding frequencies become higher), the weights $|C_k|$ decrease, indicating these high-frequency basis functions contribute less to the original function. This behavior is characteristic of a function that is more regular or smooth because it lacks high-frequency components.
The critical wavenumber $k_{cut}$ can be viewed as a turning point beyond which the magnitude of the Fourier coefficients decreases rapidly. The position of $k_{cut}$ is intrinsically linked to the parameter $n$ of the function $f(\theta) = \cos(n \sin(\theta/2) \pi)$. As $n$ increases, the function $f(\theta)$ incorporates more oscillations within the fixed interval $[0, 2\pi]$, and therefore, its Fourier transform will have significant coefficients up to a higher wavenumber. This shifts $k_{cut}$ to a larger value, implying that the function contains substantial frequency components up to this new $k_{cut}$. It is this shift that reflects the increase in $n$ and thus the increase in the number of oscillations in $f(\theta)$.
In summary, the magnitude of the Fourier coefficients $C_k$ represents the contribution of each frequency component to the original function, and a rapid decay in $|C_k|$ with respect to $k$ indicates a predominance of lower frequencies in the function's composition. Meanwhile, the dependency of $k_{cut}$ on $n$ demonstrates that more complex functions with more rapid oscillations (higher $n$) will naturally have a broader spectrum of significant frequency components before this decay commences.
The derivative $f'(t)$ is:
\begin{align*} f'(t) = -\sum_{k=1}^{\left\lfloor \frac{n}{2} \right\rfloor - 1} \left( a_k \cdot \frac{2 \pi k}{d - c} \cdot \sin\left(\frac{2 \pi k t}{d - c}\right) - b_k \cdot \frac{2 \pi k}{d - c} \cdot \cos\left(\frac{2 \pi k t}{d - c}\right) \right) \end{align*}Here's the breakdown:
This derivative captures the rate of change of the trigonometric interpolation function with respect to time $t$.
from pylab import *
import math
def get_trig_interpolation_di(dft_results, c, d, t_num=None):
if t_num is None:
n = len(dft_results)
else:
n = t_num
a0 = dft_results[0].real / n # The constant term (DC component)
an = (dft_results[1:n // 2].real * 2 / n).tolist() # The cosine coefficients
bn = (-dft_results[1:n // 2].imag * 2 / n).tolist() # The sine coefficients
# Handling the Nyquist frequency term
if n % 2 == 0:
an.append(dft_results[n // 2].real / n)
else:
an.append(dft_results[n // 2].real * 2 / n)
# Function to calculate the trigonometric interpolation at a given point t
def trig_interpolation_function_d(t):
result = 0
for k in range(1, n // 2):
result += -an[k - 1] * np.sin(2 * np.pi * k * t / (d - c)) * (2 * np.pi * k / (d - c))
result += bn[k - 1] * np.cos(2 * np.pi * k * t / (d - c)) * (2 * np.pi * k / (d - c))
if n % 2 == 0:
result += -an[-1] * np.cos(np.pi * (n // 2) * t / (d - c)) * (np.pi * (n // 2) / (d - c))
return result
return trig_interpolation_function_d
def calc_f_d(f_values, c, d, theta):
dft_result = np.fft.rfft(f_values)
trig_interpolation = get_trig_interpolation_di(dft_result, c, d, len(f_values))
derivative_values = [trig_interpolation(theta_i) for theta_i in theta]
return derivative_values
def f(theta):
return cos(4 * sin(theta/2))
def f_p(theta):
return -2 * sin(4 * sin(theta/2)) * cos(theta/2)
interval = [2**n for n in range(3,10)]
absolute_errors = []
# ---- another version of visualization ---
# for n in interval:
# theta = arange(0, 2*pi,2*pi/n)
# d_theta = 2*pi/n
# f_values = f(theta)
# fp_values = f_p(theta)
# derivative = calc_f_d(f_values, 0, 2 * np.pi, theta)
# max_error = max(abs(fp_values - derivative))
# absolute_errors.append(max_error)
# figure()
# plot(theta, derivative, label='DFT Derivative')
# title(f'Derivative of f(theta) with n={n} points')
# xlabel('Theta')
# ylabel('Derivative')
# legend()
# show()
import matplotlib.pyplot as plt
# Assuming 'interval' is a list of values for 'n'
rows = len(interval) # One row for each plot
fig, axes = plt.subplots(rows, 1, figsize=(10, 2*rows)) # Adjust figsize to make plots wider and shorter
for idx, n in enumerate(interval):
theta = np.arange(0, 2*np.pi, 2*np.pi/n)
d_theta = 2*np.pi/n
f_values = f(theta) # Make sure 'f' is defined
fp_values = f_p(theta) # Make sure 'f_p' is defined
derivative = calc_f_d(f_values, 0, 2 * np.pi, theta) # Make sure 'calc_f_d' is defined
max_error = max(abs(fp_values - derivative))
absolute_errors.append(max_error) # Make sure 'absolute_errors' is initialized
ax = axes[idx] if rows > 1 else axes # If only one plot, axes is not a list
ax.plot(theta, derivative, label='DFT Derivative')
ax.set_title(f'Derivative of f(theta) with n={n} points')
ax.set_xlabel('Theta')
ax.set_ylabel('Derivative')
ax.legend()
plt.tight_layout() # Adjust subplots to fit in the figure area
plt.show()
plot(interval, absolute_errors)
xlabel("Number of Points", fontsize=20)
ylabel("Absolute Error", fontsize=20)
title("Absolute Error")
yscale('log')
xscale('log')
Other than calculating the derivative of intepolating function, I found another approach also give the safe effect:
from pylab import *
import math
def f(theta):
return cos(4 * sin(theta/2))
def f_p(theta):
return -2 * sin(4 * sin(theta/2)) * cos(theta/2)
def dftDerivative(f_values,d_theta):
dft = rfft(f_values)
k = rfftfreq(len(f_values), d_theta)
dft_derivative = dft * (1j * 2 * pi * k)
derivative = irfft(dft_derivative, len(f_values))
return derivative
interval = [2**n for n in range(3,10)]
absolute_errors = []
import matplotlib.pyplot as plt
# Assuming 'interval' is a list of values for 'n'
rows = len(interval) # One row for each plot
fig, axes = plt.subplots(rows, 1, figsize=(10, 2*rows)) # Adjust figsize to make plots wider and shorter
for idx, n in enumerate(interval):
theta = np.arange(0, 2*np.pi, 2*np.pi/n)
d_theta = 2*np.pi/n
f_values = f(theta) # Make sure 'f' is defined
fp_values = f_p(theta) # Make sure 'f_p' is defined
derivative = calc_f_d(f_values, 0, 2 * np.pi, theta) # Make sure 'calc_f_d' is defined
max_error = max(abs(fp_values - derivative))
absolute_errors.append(max_error) # Make sure 'absolute_errors' is initialized
ax = axes[idx] if rows > 1 else axes # If only one plot, axes is not a list
ax.plot(theta, derivative, label='DFT Derivative')
ax.set_title(f'Derivative of f(theta) with n={n} points')
ax.set_xlabel('Theta')
ax.set_ylabel('Derivative')
ax.legend()
plt.tight_layout() # Adjust subplots to fit in the figure area
plt.show()
plot(interval, absolute_errors)
xlabel("Number of Points", fontsize=20)
ylabel("Absolute Error", fontsize=20)
title("Absolute Error")
yscale('log')
xscale('log')
from pylab import *
import math
import numpy as np
import matplotlib.pyplot as plt
# Define the function
def error(x, h):
correct = f_p(x)
approx = calc_diff_5p(g, x, h)
return abs(correct - approx)
def calc_diff_3p(f, x, h):
return (f(x + h) - f(x - h)) / (2 * h)
def calc_diff_5p(f, x, h):
return (f(x - 2*h) - 8 * f(x - h) + 8 * f(x + h) - f(x + 2 * h)) / (12 * h)
def g(theta):
return cos(4 * sin(theta/2))
def f_p(theta):
return -2 * sin(4 * sin(theta/2)) * cos(theta/2)
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
# Define the function
def f(x, y):
return x**2 + y**2 + 1 # Adding 1 to avoid log(0)
# Create a meshgrid of x and y values
x = np.linspace(1e-20, 2*pi, 512) # Start from 0.1 to avoid log(0)
y = np.linspace(1e-20, 1, 10000)
X, Y = np.meshgrid(x, y)
# Evaluate the function at each point
Z = error(X, Y)
# Create the heatmap with logarithmic scale for Z values
plt.figure(figsize=(8, 6))
heatmap = plt.contourf(X, Y, Z, 50, cmap='inferno', norm=LogNorm())
plt.colorbar(heatmap)
plt.xlabel('x')
plt.ylabel('h')
plt.title('Log-Scale Heatmap of error of 5 points')
# Set both axes to logarithmic scale
plt.gca().set_xscale('log')
plt.gca().set_yscale('log')
plt.show()
import plotly.graph_objects as go
import numpy as np
Z_log = np.where(Z > 0, np.log10(Z), 0)
# Assuming the error, calc_diff_5p, and g functions are defined as above
# Create the meshgrid of x and y values
x = np.linspace(1e-20, 2*np.pi, 512)
y = np.linspace(1e-20, 1, 1000)
X, Y = np.meshgrid(x, y)
# Compute Z values using your error function
Z = error(X, Y)
# Create the contour plot
fig = go.Figure(data=
go.Contour(
z=Z_log, # Use the log-transformed Z values
x=X[0], # x coordinates
y=Y[:, 0], # y coordinates
colorscale='Inferno',
colorbar=dict(
title='Log(Error)', # Update the colorbar title to reflect the log scale
tickvals=[0, 1, 2, 3], # Example tick values
ticktext=['1', '10', '100', '1000'] # Corresponding tick labels
),
contours=dict(
coloring='heatmap',
showlabels=True,
labelfont=dict(
size=12,
color='white',
),
)
)
)
# Update layout with log scale and titles
fig.update_layout(
title='Log-Scale Heatmap of error of 5 points',
xaxis_title='x',
yaxis_title='h',
xaxis_type='log',
yaxis_type='log',
width=800,
height=600
)
fig.show()